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£NJ Abstract 

(****) The aim of this paper is to propose a new numerical model to simulate 2D vesicles interacting 

^SJ with a newtonian fluid. The inextensible membrane is modeled by a chain of circular rigid 

particles which are maintained in cohesion by using two different type of forces. First, a spring 
force is imposed between neighboring particles in the chain. Second, in order to model the 

Ph bending of the membrane, each triplet of successive particles is submitted to an angular force. 

Numerical simulations of vesicles in shear flow have been run using Finite Element Method 

7— I and the FreeFem++[l] software. Exploring different ratios of inner and outer viscosities, we 

recover the well known "Tank-Treading" and "Tumbling" motions predicted by theory and ex- 
periments. Moreover, for the first time, 2D simulations of the "Vacillating-Breathing" regime 
predicted by theory in [2] and observed experimentally in [3] are done without special ingredient 
like for example thermal fluctuations used in [4] . 

Keywords: Vesicle, Penalty Method, Finite Element Method, Two-Fluid Flow, Stokes 
Ch Equations, Inelastic Contact 



X 



1. Introduction 



> 

m 

f*^**) A vesicle is a closed phospholipid membrane, separating an internal fluid from the external 

suspending medium. Studying collections of vesicles embedded in a flow is of great interest, 
^vq since such systems consist in an efficient and flexible model for much complex flows as blood 

(see [5, 6, 7] for example). Of course, the simplicity of the model cannot reproduce the com- 
plexity of the biological and biochemical activities of living cells. However, vesicles can be 
considered as a simple model to study some mechanical properties of Red Blood Cells, which 
has a great influence on their global comportment. 

From the numerical point of view, vesicle simulations have been carried out using several 
numerical methods. We can mention the molecular dynamics models [8], the boundary integral 
methods in unbounded domains [9, 10, 11] and recently in confined domains [12, 13, 14]. There 
is also numerical studies of vesicle dynamics based on finite difference methods coupled with 
level set techniques (see [15, 16] for example). 
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To the best our knowledge, finite element methods are seldom used to simulate vesicle dy- 
namics. However, we can cite [17] for a contact algorithm coupled with finite element and these 
recent works [18, 19, 20] for coupling Level Set with Finite Element. 

Apart from level set techniques, we can also mention the Phase Field method in the class 
of capture interface methods used to simulate vesicles. See for example [21, 22] or [23] for a 
numerical analysis of this method. 

We can also refer to the numerical method developed in [24] for Red Blood Cells simulations 
and where the main ingredient is to introduce two energies for elasticity and bending in order 
to model the RBC membrane. The whole problem, made of fluid and membrane interaction, 
is solved using a moving particle semi-implicit method. In the same context, we find in [4] a 
theoretical and numerical study of 2D vesicles under shear. The numerical method is based on a 
mesoscale simulation technique in which the membrane is modeled as a set of monomers con- 
nected thanks to a spring potential. The bending of the membrane is also taken into account using 
a bending potential. Theses potentials are similar to the energies used in [24]. The latter work 
is particularly interesting because it covers for the first time the Vacillating-Breathing motion in 
2D. This is done by adding a thermal fluctuation to the membrane model. 

In this paper, we study the dynamic of a 2D vesicle in shear flow by introducing a new model. 
This method consists in modelling the vesicle's membrane by a "necklace" of circular rigid 
particles which are maintained in cohesion by using a spring force imposed between neighboring 
particles in the chain. Regarding the bending of the membrane, it is modeled by imposing an 
angular force on each triplet of successive particles. These two forces are derived form energies 
which are somewhat similar to those used in [4] and [24]. 

The whole problem is solved using Finite Element Method and the FreeFem++[l] software. 
From a numerical point of view, the rigidity of each particle of the the chain is imposed using 
a penalty method, which physically consists in making the viscosity go to infinity in the rigid 
domain (see [25, 26]). In order to avoid particles overlapping, we generalize the inelastic contact 
model proposed in [26]. Finally, a constraint is added, in order to impose the vesicle to have a 
constant area. 

This paper is organized as follow : First, we introduce our model for the vesicle's membrane 
and its coupling with Stokes equations. Second, we describe the numerical method and the 
way to impose the different constraints (rigid-body motion, constraint of non overlapping and 
constant area of the vesicle). Third, we introduce some relevant physical parameters and we 
finish by some numerical results to validate our model. This validation is done by using the 
well known behaviors of vesicles : equilibrium shapes, "Tank-Treading" (TT) and "Tumbling" 
(TB) motions. Finally, we present the transition between TT and TB regimes and we recover the 
"Vacillating-Breathing" (VB) (also called "Swinging" (SW)) motion in 2D. 

2. Modeling 

Let's recall that a vesicle is a system of two fluids separated by a membrane. In this work, 
we are in two dimensions and we denote by Q a bounded domain in R 2 that represents the whole 
system (internal and external fluids and the membrane). 

We model the vesicle membrane by a "necklace" of N rigid circular particles (See figure 1). 
We denote by CB;);=i..JV these particles, x,- their centers and r their common radius. The domain 
Q. is devided into three sub-domains: B = U^jfi, is the whole rigid domain (rigid by parts), and 
£2,„ and Q ouf correspond to the internal fluid and the external suspending medium respectively. 
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These two last domains are supposed to be filled with fluids governed by Stokes equations with 
respective viscosities fi in and fi m „. 



Bi(Xi,r) 




■Oil! 



Figure 1: Necklace model of the vesicle. 



In order to model properly the membrane, let's recall the main properties of vesicles 

• Bending Energy : it costs some energy to bend the membrane. This energy has to be 
taken into account in the model 

• Local Inextensibility : the membrane has to keep the same area (perimeter in 2D) during 
the simulation time 

• Constant Volume : the vesicle membrane is often considered to be impermeable. More- 
over, inner and outer fluids are incompressibles so the vesicle volume (area in 2D) must 
stay constant. 

To satisfy the first two points, we introduce a couple of energies : 

A bending energy which consists in submitting each triplet of rigid (successive) particles to an 
angular force that tends to align them (see figure 2). Therefore, ideally (when the vesicle 
is not deflated) the shape of the membrane is circular (spherical in 3D). This energy may 
be written as 



where k Hi is the spring constant acting on (Bj-i, B,, Bj+i) triplet, e, is the unit vector between 
particles B,_i and B,- and e !+ i is the unit vector between particles B, and Bj+i. These vectors 
are given by (see figure 2) 



In the following we assume that all - angular - springs have the same constant k a . 

A stretching energy whose role is to keep sticking each pair of successive particles through an 
elastic force. We can simply write this energy in this form 




(1) 




(2) 




(3) 
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where k rp , stands for the spring constant between particles (B h B i+1 ) (cf. figure 2), £$. for 
its free-length and { t for the length of the i ,h spring. As in the case of bending energy, we 
assume that all springs have the same constant k rp . = k rp and free-length (q. - {q. 

The Constant Area is insured by imposing a constraint which will be described in section 3.2. 




Figure 2: Notations. 

The forces derived from these energies will be plugged in the right hand side of fluids model. 
This gives the major advantage to decouple entirely fluid solver and membrane model. Thus, the 
unknowns of the problem are the velocity and pressure fields in Q \ B, respectively denoted by 
u = («i,M2) and p, together with the translational and angular velocities of the rigid particles, 
denoted by V = (V u ■ ■ ■ , V N ) s R 2N and oj = (oj u . . . , oj n ) e R N . 

Most of known behaviors of vesicles under flow take place at very low Reynolds number so 
it is natural to assume that, at each time f, the fluid flow in Q \ B — Q \ B(t) is governed by Stokes 
equations with given Dirichlet boundary conditions uq on <9£X 

—/iAu + Vp = in Q \ B, 

Vk = inQ\B, (4) 
u — uq on dQ. 

where p is equal to p m in Q,„ and \x out in Q ou ,. In addition, we impose a no-slip boundary 
condition on dB: 

u(t, x) = Vi(t) + w,(?)0 - Xi(t)) 1 - on dBiff), Vi, (5) 
where x 1 - = (x I ,x 2 ) J_ is the vector (-x 2 ,xi). Finally, these equations are coupled by writing the 
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equilibrium of forces and torques on each B,: 

f B (x - XiY ■ ft - J aB (x - Xi y -an = 0, Vi, 

where fi is the volume force exerted on B, (fi = (f. st + f: )/|B;|), n the outer normal to Q \ B 
and ff* end are the forces derived from stretching and bending energies respectively. Finally, 
recall that cr stands for the stress tensor which can be written as <x = 2pD{u) - pi where we 
denoted by D(m) the deformation tensor (D(m) = v " + < v "> ) 

A first difficulty to solve this problem is the dependence in time of each part of the domain 
(£2,„, Q ou , and B,) which imposes a priori to use time dependent conform meshes. A second one 
is the explicit coupling between the membrane model and the hydrodynamic forces. We will see 
in section 3 that we can deal with these difficulties using a Penalty method. 



3. Numerical methods 



3.1. A first fluid/particle algorithm 

Instead of solving problem (4)-(5)-(6) in it's original form, we follow the method proposed in 
[25, 26]. More precisely, we use a penalty method to impose the rigid motion on B. This allows 
us to rewrite the problem initialy defined on a moving domain Q \ B(t) into a global problem 
defined on £1 Thus, the mesh of Q that we consider can be fixed once for all and does not change 
in time. Moreover, this allows us to handle the fluids/rigid particles interaction implicitly rather 
than discretizing problem (6). To do so, following a Fictitious Domain approach, we start by 
introducing the space Kg, made of functions defined on the whole domain Q, and in which we 
will search for a velocity field 

fj = (ve ffo(Q), V/ 3(V h u)i) e M 2 x R, v = V,- + co^x - Xj) 1 - a.e. in B,} . (7) 

This space can be rewritten as follow (see [25]) 

K B = [v e H l {Q), D(v) = a.e. in b} . (8) 

We see clearly in this latter expression that elements of K B represent velocity fields in D. which 
are subject to rigid body constraint in B. Moreover, it can be shown that problem (4)-(5)-(6) is 
equivalent to the following variational formulation in the whole domain Q. (see [25] for more 
details) : 

Find u in Kb and p in Lq(Q) such that 



CP) 



2p \ D(m) : D(fi) - [ pV ■ u = ( f • u, VS e Kg, 
Jn Jn Jn 



(9) 



qV-u = 0, Sqe Lq(Q), 



where L^iQ) stands for the set of functions in L 2 (Q) with zero-mean and / = ^/jlg. (1^. 

!=1 

denotes the characteristic function of B,). 
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As it is difficult to construct finite element functions in K B , we choose to use a penalty method 
to approximate the constraint problem (P) by a sequence CP £ ) of unconstrained problems. It 
reads 



Find u E in //,,(Q) and p e in Lq(Q) such that 



CP E ) 



2p f D( M £ ) : D(«) + - f D(m £ ) : D(fi) - 



Jn Jn 



/•fi, Vfi €//(',(□), 



J <?V • M £ = 0, e Lg(O). 



^ Jn 



(10) 



One can prove (see [25]) that, when s goes to zero, u B tends to the solution u of (P). Note that this 
method can also be described by considering the rigid domain as a fluid with infinite viscosity. 

Finally, let's the time step St be given and suppose that the time interval [0, T] is discretized 
using M + 1 points t° = 0, . . . , f, . . . , t M = T with f" +1 - f = St. 

For each «, if the rigid domain B" and the external force f at time f are known, we propose 
the algorithm 1 to compute the new position of B n+l . Note that u", p", V", B n+l and x" +1 



Algorithm 1 First Algorithm 
l: Compute (u",p n ) solution of (P E ) with B - B" and / = f, 

1 r 

2: Compute the corresponding velocities of the particles: V" = — — I u", 



3: Compute the center of B" +1 , using the identity x" +x = x" + StV". 



depend on the penalty parameter s. However, the superscript s have been skipped for the sake of 
readability. 

3.2. Extra constraints and final algorithm 

Two problems may occur when we solve the vesicle problem by using the previous algo- 
rithm 1. First, nothing is done to avoid overlaps between the particles and to keep the neigh- 
bouring particles of the membrane sticking together. Then, although the velocity is supposed to 
be divergence-free, numerical approximation can make the area of the vesicle change during the 
computation and one has to pay attention to this point when running long-time simulation. 

In order to ensure, at each time step, that particles do not overlap and that two neighbouring 
particles are stuck, we add a projection step to the previous algorithm. We project the velocities 
of the particles, computed from the penalized problem (P s ) onto a set of admissible velocities. 

To choose this set, we generalize the method proposed in [27] to deal with inelastic contacts. 
More precisely, let D, ; (x") be the distance between particles i and j at time n and Gjj(x") its 
gradient. We first linearize the distance Djj(x" +l ) as follow 




D u {x n+X ) = Dij(x" + StV") ^ Dijix") + StGijix") ■ V", 



(11) 



then, we chose this discrete set of admissible velocities at time n: 




Dijix") + StG u (x") ■ V > V/ < j 

Di J+l (x n ) + StG u+l (x") ■ V < V/ e [I...N -I] 

D NA ix") + StG N Aix") ■ V <0 



(12) 



6 



A second projection step is added in order to deal with the constant area constraint. To do so, 
let's assume that the area of the vesicle is approximated by the area of the polygon defined by the 
center of the particles. Then, we compute the area variation dA induced by the displacement dx, 
of particle i. We denote by nf and nj the two outer unit normals to the polygon at point X; (see 
figure 3) and we define n, = - Xi\n* + \Xi-\ - Xj\nJ) (note that indices has to be adapted for 

i = 1 and i = N). Then, the considered area variation is given by 




dA = dA i+ i + dAi_\ 

\x M - x\dxi ■ n\ \xi-\ - Xi\d,Xi ■ nT 



2 

— rii ■ dxi. 

Obviously, in the case where all particles move we obtain 

N 



dA = rii ■ dxi, 



and the constraint at time n is written 

N 



A = A n+l =A"+dA=A n + dtj^rii ■ Vf. 



!=1 

Finally, the space of admissible velocities for the area constraint at time n is given by 



= e R 2N , A = A" + St J] m ■ Vi J 



and the global algorithm is given by algorithm 2, where 11^ denotes the projection onto K and 
is performed using a Uzawa algorithm. 
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Algorithm 2 Final Algorithm 
l: Compute (u",p") solution of (!P e ) with B = B" and f = /", 

1 r 

2: Compute the corresponding velocities of the particles: V" = I u , 

3: Deal with contacts by projecting V" on /f" : V" = n K »V", 

4: Deal with the area constraint by projecting V" on K" : V" = Tl^ V", 

5: Compute the new position of B" +l : x' t ' +l — x" + 6tV", 



4. Dimensionless numbers 

To validate our model by comparing its results with those from the literature, we need to 
define some relevant dimensionless parameters which control the problem. 

The reduced Area a is defined as the ratio between the vesicle's area and the area of a circle 
having the same perimeter. 

a=4>' (13) 

where A is the area of the vesicle and Rq its effective radius given by £- (P being the 
perimeter of the vesicle). 

The reduced area a measure the amount of deflation of the vesicle. We will see in section 
5 that this parameter controls the equilibrium shape of the vesicle. 

The Viscosity Contrast A is defined as the ratio between the viscosities of the inner and the 
outer fluids. It reads 

A=^. (14) 

Pout 

This parameter monitors the dynamic of the vesicle (see section 5). 
The confinement t is given by 

= 

T / 

where / is the height of the chanel in which the vesicle evolves. 
The Capilary number C a (see Appendix A) which can be defined as 

C " - ~2k7~ (16) 

where y stands for the shear rate and r for the radius of the particles in our model. 

The Capillary number can be seen as the ratio between the characteristic time of the shear 
and a characteristic time related to the bending force. This number measures the amount of 
deformation of the vesicle. Indeed, for high values of C a , the vesicle is more deformable 
due to the hydrodynamic forces, while for low C a , the bending energy requires a higher 
cost to change the curvature of the vesicle. 



(15) 
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5. Numerical Results 



We present in this section some numerical results obtained by using algorithm 2. The penal- 
ized Stokes problem (!P £ ) is discretized using first order Finite Element where the discrete space 
for the pressure is based on Lagrange polynomials P l and the discrete space for the velocity field 
is based on P l bubble (which is an enrichment of P l to satisfy the discrete inf-sup condition). 

For each of the following simulations, the time step is 6t = 5.1CT 3 , the reference length of 
the springs is €q = 2r where r — 1 .5 is the radius of the rigid particles, the penalty parameter is 
s = 5.1CT 3 and the shear rate is y — 1. 



The other parameters are, for the different sections: 



Section 


N 


k a 




/ 


L 


Pin 


P-out 


Equilibrium shapes 


42 


200 


0.25 


150 


150 


1 


1 


TT 


50 


600 


0.5 


242 


300 


1 


1 


TB 


50 


600 


0.5 


242 


300 


1 




VB 


50 


600 


0.5 


242 


300 


1 





5.1. Equilibrium shapes 

An essential step for the validation of any vesicle model is to retrieve the equilibrium shape 
of a vesicle in a fluid at rest. More precisely, if we consider a vesicle in a shape which is not the 
one minimizing its energy, the curvature force makes the membrane move as the fluid around it. 
This dynamic stops when the vesicle reaches its equilibrium state in which the final shape of the 
vesicle depends on its reduced area. 

In our simulations the initial distribution of the rigid particles (that form the vesicle mem- 
brane) is done by placing them in an ellipse. Then, the Stokes problem is solved in a rectangular 
domain containing the vesicle with homogeneous Dirichlet boundary conditions. 

Figure 4 shows the evolution of a vesicle suspended in a fluid at rest by plotting the position 
of each rigid particle center. One can see that the vesicle reaches its equilibrium shape which is 
biconcave in the case of small reduced area. Note that this form is typical of red blood cells 
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Figure 4: Equlibrium shapes. Shapes versus t for a = 0.42. 
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It is important to check that the area and the perimeter are well conserved. During this 
simulation, the variation of the area ranges from -0.3% to 0.2% and the variation of the perimeter 
from -0.05% to 0.05% (see figure 5). 




Figure 5: Equilibrium shapes. Area and perimeter versus t for a = 0.42. 



Figure 6 shows the computed equilibrium shapes of vesicles with different reduced area. As 
it is well known, the equilibre shape goes from circular one when a — 1 to biconcave one when 
a is small enough (about 0.6). 
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Figure 6: Equilibrium shapes for different values of a. 



Finally, we compare our results to those presented in [28] (in this paper equilibrium shapes 
are computed using two numerical methods : Boundary Integral and Lattice Boltzman). Figure 
7 shows at the same times our results and those from [28] for different value of the reduced area 
parameter. One can notice that there is a good agreement between these results. Nevertheless, 
we see small differences in the cases where a is small and we believe this is due to the thickness 
of our membrane. 
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a = 0.42 



<y = 0.51 




5.2. "Tank-Treading" Regime (TT) 

We are interested in this section to the dynamic of a vesicle under linear shear flow. It is well 
known that if A, the ratio between the inner and the outer fluid viscosities, is smaller than a certain 
critical value (around 5 in unbounded domain [9]), the vesicle takes a steady angle with respect 
to the horizontal axis and its membrane starts to rotate like a chain of a tank (this motion is called 
Tank-Treading). To validate our model in this regime, we refer to the results given in [9]. So, we 
consider vesicles of different reduced area subject to linear shear flow in a rectangular domain 
and the confinement parameter t is chosen to be small (around 0.2) to minimize the effect of 
walls. In figure 8 we plot the steady angle for different reduced area and we compare our results 
with those from [9] . 
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Figure 8: Tank treading angles. Comparison between our results (in the case of small confinement r = 0.2 ) and those 
from [9] (in unbounded domain). 

Finally we plot in figure 9 the streamlines of the velocity field, its magnitude and the position 
of the vesicle at different times. In this simulation the reduced area is a — 0.85, the confinement 
parameter is r = 0.32, the viscosity contrast A = 1 and the number of rigid particles is N - 38. 

5.3. Tumbling Regime (TB) 

When the viscosity contrast A is above a certain critical value Ac, a vesicle subject to linear 
shear flow starts to rotate and tends to follow a solid rotation. This motion is called tumbling and 
is also typical of red blood cells. 

Figure 10 shows the streamlines of the velocity field, its magnitude and the position of the 
vesicle at different times. The initial time corresponds to the first time step iteration and the final 
one corresponds to the time during which the vesicle made half a turn. In this simulation we 
keep the same parameters as in the case of TT simulation except the viscosity contrast A = 20. 

If we look at the particle marked in red, we note that the membrane has also turned slightly 
on itself. The vesicle rotation is not yet that of a solid. See [20] in which the authors showed 
that increasing the viscosity ratio increases the rotation frequency of the vesicle. This frequency 
reaches a steady value which corresponds to the rotation of a solid object. In this work, we show 
in figure 11 the angles of a vesicle in tumbling regime versus time for two values of A (12 and 
20). We can see easily that the frequency rotation when A = 20 is greater than in the case A = 12. 
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velocity Magnitude velocity Magnitude 




(a) r = 5.1(r 3 (b) t = 0.5 



velocity Magnitude velocity Magnitude 




(e) / = 15.7 



Figure 9: Streamlines of the velocity field and vesicle position in Tank-Treading regime at different time, a = 0.85, 
r = 0.32,JV = 38an(U=l. 

5.4. Transition between Tank-Treading and Tumbling. The "Vacillating -Breathing" Regime 
(VB) 

In addition to the two previous regimes TT and TB, it has been shown (see [2] and [3]) that 
in a shear flow, a vesicle may exhibit another kind of motion called "Vacillating-Breathing" in 
[2] and "Swinging" in [3]. This regime occurs when the vesicle's membrane is fairly deformable 
(for a high enough capillary number C„) and precedes the TB regime when increasing the vis- 
cosity ratio A. Thus, the vesicle shape undergoes large deformations, the long main axis of the 
vesicle performs oscillations around the flow direction. More precisely, f i) when the vesicle is 
perpendicular to the shear axis, it is being deformed till it attains a quasi circular shape, and then 
(ii) when it tends to be aligned with the shear axis it is elongated like in the case of TT regime. 

In addition to theoretical and experimental results presented respectively in [2] and [3] to 
describe this motion, we can find some numerical studies highlighting the VB regime but they 
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Figure 10: Streamlines of the velocity field and vesicle position in Tumbling regime at different time, a = 0.85, r = 0.32, 
N = 38 and A = 20. 



are still infrequent. See [29] for 3D simulation using Boundary Integral Method and [4] for 
2D simulation by adding thermal fluctuation. To our knowledge, our work is the first one that 
presents 2D numerical simulation of VB regime without adding thermal fluctuation. Moreover, 
in [11] and [4] the authors argue that there is no support to the existence of VB regime in 2D 
unless the addition of high-order undulation modes. 

Figure 12 shows the streamlines of the velocity field, its magnitude and the position of the 
vesicle at different times. The initial time corresponds to the first time step iteration and the final 
one corresponds to the time during which the vesicle made half a turn. In this simulation we take 
A = 7.5 as viscosity ratio and we keep unchanged all other parameters as in the cases of TT and 
TB simulations (the capillary number is then C a 7.68). 

In order to point out the fact that VB is an intermediate regime between TT and TB, we 
present in figure 14 the vesicle angle versus time for different viscosity ratios. We see clearly 
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Figure 1 1 : Angles versus time in tumbling regime, a = 0.85 and t = 0.2. 



that in this case, we have TT regime up to the critical value A — 5.5 for which the steady angle is 
close to zero. Then, for A = 6.5 we have a VB motion. Finally TB regime starts from A = 8. 

For the same purpose, we plot in figure 14 the steady TT angle versus the viscosity ratio for 
two values of the reduced area. One can see that the presence of the VB regime suppresses the 
square root singularity at the tumbling threshold as it is explained in [11] 

Finally, we plot in figure 15 the area and perimeter variations in the case of VB regime. Note 
that we have the same results in the cases of TT and TB regimes but we choose to show only the 
VB case since it is the less favorable regime because the vesicle undergoes large deformations. 



6. Conclusion 

We introduced in this work a new model to study the dynamics of 2D vesicles under shear. 

First, we described our model and pointed out the way to solve it using finite elements cou- 
pled with a penalty method. Mainly, this model is based on Stokes equations for the inner and 
outer fluids and on a necklace of rigid particles for the vesicle membrane. The physical properties 
of the vesicle are taken into account by a couple of forces acting on the rigid particles to model 
the inextensibility and the bending of the membrane. The area and perimeter conservation are 
obtained with a good precision by projecting the velocity of the particles onto a set of admissible 
fields satisfying these constraints. These projections are done using Uzawa algorithm. 

One of the advantages of our model is its link with fluid/rigid particles problems. One can 
see that it is slightly easy to extend any fluid/rigid particles solver to deal with several vesicles. 
However, it is not easy to extend it to 3D vesicles unless to authorize the overlapping of the rigid 
particles which introduces technical difficulties. 

Second, we have validated our model by retrieving the equilibrium shapes of vesicles in a 
fluid at rest and by recovering the TT and TB regimes. 

15 



Finally, we pointed out the transitional regime (when increasing the viscosity ratio) which is 
called Vacillating-Breathing or Swinging. This particular motion has been predicted theoretically 
in [2] and observed experimentally in [3]. It has been carried out numerically for 3D vesicles in 
[29] and for 2D vesicles in [4] by adding thermal fluctuations to the membrane model. To the 
best of our knowledge, our results are the first ones recovering numerically the VB regime in 2D 
without adding any special ingredient. 
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Appendix A. Numerical bending modulus and capillary number 

In this section we derive the numerical capillary number C% which corresponds to a vesicle's 
membrane represented by N rigid particles. Let us recall the bending energy E/, we used in our 
model 

E b = J^k a (e r e i+l + l). (A.l) 



We assume that the vesicle has a circular shape with radius R. For large N, namely for small 

qN _ In 
( — N ' 



6? = Tj, see figure 2 for notations, we have 



E N h = fe fl 2(cos(2 J 8f) + l)=^2(cos(^-0f) + l), 

i i 

= -& o £(cos(0f)-l). 



Using the fact that dff w where r N — ^ is the radius of the particles, one gets for large 



N : 



i 

K2r N 



2\ 



2 R 2 



1 k a 2r A 



C 1 



2 ^ ' R 2 2 

i 

From this, together with the fact that the continuous bending energy is given by 



where B is the bending modulus, we define the discrete bending modulus 

N _ k a 2nR 



B N = k a 2r n 



N 



The continuous capillary number is defined by (see [24]) : 
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B 



where is the viscosity of the outer fluid and I is the caracteristic length of the vesicle. If one 
takes I = Rwe obtain 

_ floury 

We can then define a discrete capillary number : 

c n _ Vo UI R 3 y _ Ho U iR 3 y 
a ~ B N ~ k a 2r N ' 

Note that if we let k a depend on N and if we chose 
we obtain, for all N 

B N - B and C N a = C a . 
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Figure 12: Streamlines of the velocity field and vesicle position in Vacillating-Breathing regime at different time, a = 
0.85, r = 0.32, N = 38 and A = 7.5. 
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Figure 13: Angles versus time for different viscosity ratios, a = 0.85, r = 0.2. 
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Figure 14: TT Angles versus viscosity ration. 
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Figure 15: Vacillating-Breathing. Area and perimeter versus t for a = 0.85. 
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